!	** Simplex Simulated Annealing module **
!	Alpha arrays are not stored locally here - create them in the main code and pass them through
!	Implement the cost() function as a real*8 external subroutine.

	module simplex
	implicit none
	! Fit variables
	integer :: simplex_nalpha = 0, simplex_npasses = 10, simplex_nloops = 100, simplex_ncycles = 1000
	real*8, allocatable :: simplex_alpha(:)
	contains

	subroutine simplex_createarrays()
	implicit none
	if (simplex_nalpha.lt.1) stop "Must set a useful simplex_nalpha value before calling simplex_initarrays."
	allocate(simplex_alpha(simplex_nalpha))
	end subroutine 

	subroutine simplex_freearrays()
	implicit none
	deallocate(simplex_alpha)
	end subroutine 

	! Simplex Simulated Annealing following original Nelder-Mead simplex algorithm
	real*8 function simplex_minimise()
	use random
	implicit none
	integer :: a, n, m, r, vbest, vworst, vnextworst
! 	real*8, intent(inout) :: alpha(nalpha)
	real*8 :: vertex(simplex_nalpha+1,simplex_nalpha), coldcosts(simplex_nalpha+1)
	real*8 :: vert_r(simplex_nalpha), vert_e(simplex_nalpha), vert_c(simplex_nalpha), xbar(simplex_nalpha)
	real*8 :: cost_r, cost_e, cost_c, srho, schi, ssigma, sgamma, t, cost_best, cost_worst, cost_nextworst, tempcost
	real*8 :: lscale = 0.5d0, ctol = 1.0d-6, lastbest = 0.0d0
	logical :: converged
	real*8, external :: cost
	! Simplex minimisation using original Nelder-Mead Simplex move coefficients
	srho = 1.0d0
	schi = 2.0d0
	sgamma = 0.5d0
	ssigma = 0.5d0
	!
	! Begin Master Simplex Loop
	!
	t = cost(simplex_nalpha,simplex_alpha) * 10.0
	do r = 0,simplex_nloops
	  ! Reduce temperature if not first cycle
	  if (r.ne.0) t = t * 0.5
	  write(6,"(a,f)") "           Vertex temperature for this pass is ",t
	  ! Set up primary vertices
	  do a=1,simplex_nalpha
	    vertex(1,a) = simplex_alpha(a)
	  end do
	  ! Set up remaining vertices
	  do n=2,simplex_nalpha+1
	    do a=1,simplex_nalpha
	      ! Randomise each point in this vertex a little bit 
	      vertex(n,a) = simplex_alpha(a) * ((1.0d0 - lscale) + 2.0d0 * lscale * ran())
	    end do
	  end do
	  ! Evaluate cost function at each vertex
	  do n=1,simplex_nalpha+1
	    coldcosts(n) = cost(simplex_nalpha,vertex(n,:))
	    !write(6,"('Vertex ',i3,' value = ',e12.6,20f7.3)") n,coldcosts(n),(vertex(n,a),a=1,nalpha)
	  end do
10	  format (a,i3,a,e12.6,a,e12.6)
	  write(6,10) "Pass ",r," : Vertices have 'cool' cost range of ",minval(coldcosts)," to ",maxval(coldcosts)
	  !
	  ! Begin Simplex Minimisation
	  !
	  do n=1,simplex_ncycles

	    ! Sort initial vertices so lowest cost is in vertex(1,:)
	    call sort_simplex(simplex_nalpha,vertex,coldcosts)
	    ! 1 - Find best, worst, and next-worst points (inc. random T fluctuations)
	    vbest = 1
	    cost_best = coldcosts(1)
	    vworst = 2
	    cost_worst = coldcosts(2)
	    if (simplex_nalpha.eq.1) then
	      vnextworst = 2
	      cost_nextworst = coldcosts(2)
	    else
	      vnextworst = 3
	      cost_nextworst = coldcosts(3)
	    end if
	    do a=1,simplex_nalpha+1
	      ! Get cost of vertex 'a' and do Boltzmann fluctuation
	      tempcost = coldcosts(a) - t * dlog(ran())
	      if (tempcost.lt.cost_best) then
		vbest = a
		cost_best = tempcost
	      else if (tempcost.gt.cost_worst) then
		vworst = a
		cost_worst = tempcost
	      else if (tempcost.gt.cost_nextworst) then
		vnextworst = a
		cost_nextworst = tempcost
	      end if
	    end do
	!write(0,*) vbest,vworst,vnextworst
	!write(0,*) cost_best, cost_worst, cost_nextworst
	    !write(u,"('Best : ',30f5.2)") (vertex(1,a)-10.0d0,a=1,nalpha)
	    ! ... Calculate xbar for use in simplex move routines
	    xbar = sum(vertex(a,:)) / real(simplex_nalpha)
	    ! 2 - Calculate reflection vertex
	    call simplex_reflect(simplex_nalpha,vertex,xbar,vert_r,srho)
	    cost_r = cost(simplex_nalpha,vert_r)
	    !write(0,"('R ',e10.2,40f6.2)") cost_r,vert_r
	    ! ... If this point is lower in cost than the hext-highest vertex, but higher than the best point, accept it and terminate iteration
	    if ((cost_best.le.cost_r).and.(cost_r.lt.cost_nextworst)) then
	      !write(0,*) "Reflection is better than n'th vertex, but worse than 1st, so accepting reflection."
	      vertex(simplex_nalpha+1,:) = vert_r(:)
	      coldcosts(simplex_nalpha+1) = cost_r
	      ! Check convergence, otherwise cycle to next iteration
	      if (converged(simplex_nalpha,coldcosts,ctol,n)) exit
	      cycle
	    end if
	    ! 3 - Calculate expansion point (if reflected point is lower than the best point already found)
	    ! 4 - Calculate contraction point (if reflected point is worse than the next worst (n-1'th) point)
	    if (cost_r.lt.cost_best) then
	      call simplex_expand(simplex_nalpha,vertex,xbar,vert_e,srho,schi)
	      cost_e = cost(simplex_nalpha,vert_e)
	      !write(0,"('E ',e10.2,40f6.2)") cost_e,vert_e
	      ! ... If expanded point is better than reflected point, accept expanded point and terminate iteration.
	      ! ... If reflected point is better (or equal to) expanded point, accept reflected point and terminate iteration
	      if (cost_e.lt.cost_r) then
	        !write(0,*) "Expansion is better than reflection, so accepting expansion."
	        vertex(simplex_nalpha+1,:) = vert_e(:)
	        coldcosts(simplex_nalpha+1) = cost_e
	        ! Check convergence, otherwise cycle to next iteration
	        if (converged(simplex_nalpha,coldcosts,ctol,n)) exit
	        cycle
	      else if (cost_e.ge.cost_r) then
	        !write(0,*) "Expansion is worse than reflection, so accepting reflection."
	        vertex(simplex_nalpha+1,:) = vert_r(:)
	        coldcosts(simplex_nalpha+1) = cost_r
	        ! Check convergence, otherwise cycle to next iteration
	        if (converged(simplex_nalpha,coldcosts,ctol,n)) exit
	        cycle
	      end if
	    else if (cost_r.ge.cost_nextworst) then
	      ! ... If reflected point is bettern than worst point, but worse than next-best (n-1) point, do 'outside' contraction
	      if ((cost_r.lt.cost_worst).and.(cost_nextworst.lt.cost_r)) then
	        call simplex_contract_out(simplex_nalpha,vertex,xbar,vert_c,srho,sgamma)
	        cost_c = cost(simplex_nalpha,vert_c)
	        !write(0,"('Co ',e10.2,40f6.2)") cost_c,vert_c
	        if (cost_c.le.cost_r) then
	  	!write(0,*) "Contraction(outside) is better than reflection, so accepting contraction."
	          vertex(simplex_nalpha+1,:) = vert_c(:)
	          coldcosts(simplex_nalpha+1) = cost_c
		  ! Check convergence, otherwise cycle to next iteration
		  if (converged(simplex_nalpha,coldcosts,ctol,n)) exit
	          cycle
	        end if
	      else  ! ... do an inside contraction
	        call simplex_contract_in(simplex_nalpha,vertex,xbar,vert_c,sgamma)
	        cost_c = cost(simplex_nalpha,vert_c)
	        !write(0,"('Ci ',e10.2,40f6.2)") cost_c,vert_c
	        if (cost_c.le.coldcosts(simplex_nalpha+1)) then
	  	!write(0,*) "Contraction(inside) is better than worst point, so accepting contraction."
	          vertex(simplex_nalpha,:) = vert_c(:)
	          coldcosts(simplex_nalpha) = cost_c
		  ! Check convergence, otherwise cycle to next iteration
		  if (converged(simplex_nalpha,coldcosts,ctol,n)) exit
	          cycle
	        end if
	      end if
	    end if
	    ! 5 - Shrink (if contractions fail)
	    call simplex_shrink(simplex_nalpha,vertex,ssigma)
	    ! Check convergence
	    if (converged(simplex_nalpha,coldcosts,ctol,n)) exit

	  end do

	  ! Sort vertices so lowest cost is in vertex(1,:)
	  call sort_simplex(simplex_nalpha,vertex,coldcosts)

	  ! Recalculate cost for the best vertex and print cost per species per atom decomposition
	  cost_c = cost(simplex_nalpha,vertex(1,:))
	  write(6,"(8x,a,e12.6)") "   Best fit achieved has cost of ",cost_c
20	  format (a,i3,a,i2,a,e12.6)
21	  format (12x,a,e12.6,a)

	  ! Load best parameters into alpha array
	  simplex_alpha(:) = vertex(1,:)

	end do ! Restart loop

	! Grab current best parameters
	simplex_alpha(:) = vertex(1,:)
	coldcosts(1) = cost(simplex_nalpha,simplex_alpha(:))
	
	simplex_minimise = coldcosts(1)

	end function simplex_minimise

	logical function converged(nalpha, costs, ctol, cycleno)
	implicit none
	integer :: a, cycleno, nalpha
	real*8 :: costs(nalpha+1), vmean, serror, ctol
	! Test simplex convergence
	vmean = sum(costs) / real(nalpha+1)
	serror = 0.0d0
	do a=1,nalpha+1
	  serror = serror + (costs(a) - vmean)**2
	end do
	serror = sqrt( serror / real(nalpha+1) )
	if (serror.lt.ctol) then
	  write(6,"(a,i5)") " Simplex converged at step ",cycleno
	  converged = .true.
	else
	  converged = .false.
	end if
	end function converged

	subroutine sort_simplex(nalpha,vertex,costs)
	implicit none
	integer :: nalpha, incrmnt, i, j, incr
	real*8 :: vertex(nalpha+1,nalpha), costs(nalpha+1), tempv(nalpha), tempc
	! Shell Sort - Reformatted algorithm AS 304.8 APPL.STATIST. (1996), VOL.45, NO.3

	! Set initial increment 
        incrmnt = 1
	do while (incrmnt.le.nalpha+1)
	  incrmnt = 3 * incrmnt + 1
	end do
	!write(0,*) incrmnt

	do while (incrmnt.gt.1)
	  incrmnt = incrmnt / 3
	  i = incrmnt + 1

	  do while (i.le.nalpha+1)
            tempc = costs(i)
	    tempv(:) = vertex(i,:)
            j = i
	    do while (costs(j-incrmnt).ge.tempc)
              costs(j) = costs(j-incrmnt)
              vertex(j,:) = vertex(j-incrmnt,:)
              j = j - incrmnt
              if (j.le.incrmnt) exit
	    end do
            costs(j) = tempc
	    vertex(j,:) = tempv(:)
            i = i + 1
	  end do

	end do

	end subroutine sort_simplex

	subroutine simplex_reflect(nalpha,vertex,xbar,xnew,rho)
	implicit none
	integer :: n, nalpha
	real*8 :: vertex(nalpha+1,nalpha), xnew(nalpha), xbar(nalpha), rho
	! Reflect highest point : x(new) = (1 + rho) * y - rho * x(n+1)    where y = sum(i=1,n) x(i) / n
	! Calculate x(new)
	xnew(:) = (1.0d0 + rho) * xbar(:) - rho * vertex(nalpha+1,:)
	end subroutine simplex_reflect

	subroutine simplex_expand(nalpha,vertex,xbar,xnew,rho,chi)
	implicit none
	integer :: nalpha
	real*8 :: vertex(nalpha+1,nalpha), xnew(nalpha), xbar(nalpha), rho, chi
	! Expand around highest point : x(new) = (1 + rho * chi) * xbar - rho * chi * x(n+1)    where xbar = sum(i=1,n) x(i) / n
	xnew = (1.0d0 + rho * chi) * xbar - rho * chi * vertex(nalpha+1,:)
	end subroutine simplex_expand

	subroutine simplex_contract_out(nalpha,vertex,xbar,xnew,rho,gamma)
	implicit none
	integer :: nalpha
	real*8 :: vertex(nalpha+1,nalpha), xnew(nalpha), xbar(nalpha), rho, gamma
	! Contract simplex (new point outside current polytope) : x(new) = (1 + rho * gamma) * xbar - rho * gamma * vertex(nalpha+1)
	xnew = (1.0d0 + rho * gamma) * xbar - rho * gamma * vertex(nalpha+1,:)
	end subroutine simplex_contract_out

	subroutine simplex_contract_in(nalpha,vertex,xbar,xnew,gamma)
	implicit none
	integer :: nalpha
	real*8 :: vertex(nalpha+1,nalpha), xnew(nalpha), xbar(nalpha), gamma
	! Contract simplex (new point inside current polytope) : x(new) = (1 - gamma) * xbar + gamma * vertex(nalpha+1)
	xnew = (1.0d0 - gamma) * xbar + gamma * vertex(nalpha+1,:)
	end subroutine simplex_contract_in

	subroutine simplex_shrink(nalpha,vertex,sigma)
	implicit none
	integer :: n,nalpha
	real*8 :: vertex(nalpha+1,nalpha), xnew(nalpha), sigma
	! Shrink vertices of current simplex, leaving only the first (x(1)) as-is: x(n) = x(1) + sigma(x(n) - x(1)),   n=2,nalpha+1
	do n=2,nalpha+1
	  vertex(n,:) = sigma * ( vertex(n,:) - vertex(1,:) ) + vertex(1,:)
	end do
	end subroutine simplex_shrink

	end module simplex
